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ABSTRACT 

Star-forming galaxies which are too faint to be detected individually produce intensity 
fluctuations in the cosmic background light. This contribution needs to be taken into 
■ account as a foreground when using the primordial signal to constrain cosmological pa- 

rameters. The extragalactic fluctuations are also interesting in their own right as they 
depend on the star formation history of the Universe and the way in which this con- 
nects with the formation of cosmic structure. We present a new framework which allows 
us to predict the occupation of dark matter haloes by star-forming galaxies and uses 
this information, in conjunction with an N-body simulation of structure formation, to 
predict the power spectrum of intensity fluctuations in the infrared background. We 
compute the emission from galaxies at far-infrared, millimetre and radio wavelengths. 
Our method gives accurate predictions for the clustering of galaxies both for the one 
1 halo and two halo terms. We illustrate our new framework using a previously pub- 

lished model which reproduces the number counts and redshift distribution of galaxies 
selected by their emission at 850 fim. Without adjusting any of the model parameters, 
the predictions show encouraging agreement at high frequencies and on small angular 
scales with recent estimates of the extragalactic fluctuations in the background made 
from early data analysed by the Planck Collaboration. There are, however, substantial 
discrepancies between the model predictions and observations on large angular scales 
and at low frequencies, which illustrates the usefulness of the intensity fluctuations as 
a constraint on galaxy formation models. 
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1 INTRODUCTION to be removed statistically in order to get to the primordial 

CBL signal, or as an interesting quantity in its own right as 
The cosmic background light (CBL) is a rich source of lnlor- , „ , , r , , . r , , . 

, , , , a challenge to models ot the clustering ol galaxies and their 

mation about the conditions m the early universe and the . . . „ . , , , ,. 

emission m the mlra-red, millimetre and radio ranges ot the 
subsequent growth ol galaxies and ol structure m the dark , 

electromagnetic spectrum. 

matter. Accurate measurements ol the power spectrum ol 

temperature anisotropies in the primordial component have A sma11 fraction, around 10-20%, of the extragalactic 

led to tight constraints being placed on the values of the background light has been resolved into galaxies at far- 
basic cosmological parameters (e.g. Komatsu et al. 2011). infrared and millimetre wavelengths (Bethermin et al. 2010; 
Other contributions to the CBL include Galactic cirrus, the ° liver et aL 2011 )- Fluctuations in the intensity of the unre- 
thermal and kinetic Sunyaev-Zel'dovich effects and extra- solved CBL have been discovered recently (Grossan & Smoot 
galactic sources such as star-forming galaxies (SFGs) and 2007 ' Lagache et al. 2007; Viero et al. 2009; Hall et al. 2010; 
active radio galaxies (see e.g. Figure 2 of Dunkley et al. the Planck Collaboration 2011b; Penin et al. 2012a). These 
2011). Correlated fluctuations in the CBL due to SFGs de- fluctuations have two sources: the shot noise arising from 
pend on the number of sources and their clustering, which sampling a discrete number of unresolved galaxies within 
in turn is sensitive to the variation in the efficiency of star the telescope beam and the intrinsic clustering of the galax- 
formation with dark matter halo mass. The extragalactic ies - In an earl y an alysis of six regions of low Galactic extinc- 
contribution to the CBL may be viewed either as a nuisance tion coverin g 140 square degrees, the Planck Collaboration 

(2011b) have cleaned the temperature anisotropy maps to 
leave only the extragalactic fluctuations in the CBL. This is 
* hansikk@unimelb.edu.au done by using the lowest frequency Planck map to remove 
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the cosmological signal (we shall see later that fluctuations 
due to the clustering of extragalactic sources are negligible 
at low frequencies) and exploiting neutral hydrogen obser- 
vations as a tracer of dust to further reduce the contribution 
from Galactic emission. (For an overview of the Planck mis- 
sion, see the Planck Collaboration 2011a.) 

A variety of models have been developed to interpret 
the measured power spectrum of fluctuations in the inten- 
sity of the CBL. These models have predominantly used 
empirical spectral energy distributions (e.g. Lagache et al. 
2003, 2007). Simple analytic models have been assumed for 
the clustering of galaxies such as a linear bias factor rela- 
tive to the clustering of the dark matter (e.g. Knox et al. 
2001; Fernandez-Conde et al. 2008; Hall et al. 2010) or the 
halo occupation distribution (HOD) formalism (Amblard & 
Cooray 2007; Viero et al. 2009; Shang et al. 2011; Penin et al. 
2012b). Righi et al. (2008) presented a calculation based on 
the mergers of dark matter haloes and a simple dust evolu- 
tion model. Seghal et al. (2010) assumed that the number of 
infra-red sources scaled with halo mass to compute the two 
halo clustering, without calculating a one halo contribution. 
Negrello et al. (2007) combined the model of Granato et al. 
(2004) for the evolution of the spheroid population with a 
phenomenological model for the evolution of starbursts, nor- 
mal late- type spirals and radio galaxies. These authors as- 
sumed a linear bias to model the clustering of galaxies. The 
Planck Collaboration (2011b) have used their measurements 
of the CBL fluctuations to rule out a linear bias model that is 
constrained to match the observed number counts of galax- 
ies, arguing that accurate small-scale clustering predictions 
are critical to match the observations. 

In this paper we present a new approach for comput- 
ing the contribution of SFGs to the intensity fluctuations 
in the CBL, with two important improvements over previ- 
ous theoretical models. First, we make an ab initio calcula- 
tion of the spectral energy distributions of a large sample of 
galaxies, using a self-consistent treatment of the extinction 
of starlight by dust and the reprocessing of the absorbed 
energy to longer wavelengths. We also compute the radio 
emission from star-forming galaxies (Condon 1992; Bressan, 
Silva & Granato 2002). Second, we combine the predictions 
for the properties of the galaxy population with a high res- 
olution, large volume N-body simulation of the clustering of 
matter in the Universe. This allows us to accurately model 
the clustering of galaxies over a wide range of pair separa- 
tions, including the highly nonlinear regime corresponding 
to scales within an individual dark matter halo. Empirical 
models suffer from the obvious drawback that the bulk of 
the galaxies responsible for the extragalactic background 
have not yet been observed, which makes the calibration 
of this class of model uncertain. Also, without the context 
of a model for structure formation, any assumptions about 
the clustering of the galaxies are decoupled from their abun- 
dance (e.g. as in the calculation by Xia et al. 2012). 

The first step in our calculation is to generate predic- 
tions for the star formation histories of a representative sam- 
ple of galaxies, using the semi-analytical galaxy formation 
code, GALFORM (Cole et al. 2000; Baugh et al. 2005). The 
star formation and metal enrichment histories, along with 
the size of the disc and bulge components are input to the 
spectro-photometric code GRASIL (Silva et al. 1998). GRASIL 
is used to compute the spectral energy distribution of each 



galaxy, using a radiative transfer calculation in a two-phase 
dust medium (Granato et al. 2000). We describe how the 
hybrid GALFORM plus GRASIL code is implanted into a large- 
volume, high resolution N-body simulation of the clustering 
of the dark matter to add information about the spatial dis- 
tribution of galaxies. This allows an accurate treatment of 
the one-halo term and of the nonlinear component of the 
two-halo term. 

The content of the paper is as follows. In Section 2 we 
give a brief overview of the GALFORM and GRASIL models be- 
fore explaining how we populate an N-body simulation with 
model galaxies. In Section 2.3 we set out the calculation of 
the angular correlation function of flux. The results of the 
paper are presented in Section 3 and the summary in Sec- 
tion 4. The appendix discusses the sensitivity of the model 
predictions to the finite resolution of the N-body simulation. 



2 THEORETICAL BACKGROUND 

In this section we introduce the galaxy formation model used 
and outline the theoretical concepts needed in the paper. We 
give a brief overview of the semi-analytical galaxy formation 
model and explain how the emission from dust and at radio 
wavelengths is computed in Section 2.1. The implementa- 
tion of this model in an N-body simulation is described in 
Section 2.2. In Section 2.3 we set out the equations describ- 
ing how the clustering of intensity fluctuations due to ex- 
tragalactic sources is computed from the predictions of the 
galaxy formation model. 

2.1 The hybrid galaxy formation model 

We use a hybrid model consisting of the GALFORM semi- 
analytical galaxy formation code (Sec 2.1.1) and the GRASIL 
spectrophotometric code (Sec 2.1.2). 

2.1.1 The GALFORM galaxy formation model 

The formation and evolution of galaxies within the ACDM 
cosmology is predicted using the semi-analytical model 
GALFORM (Cole et al. 2000). The main processes modelled 
include: (1) the formation of dark matter haloes by mergers 
and the accretion of smaller haloes, (2) the growth of galac- 
tic discs following the shock-heating and radiative cooling of 
gas inside dark matter haloes, (3) star formation in galactic 
discs, (4) the reheating and ejection of gas by supernovae, (5) 
the prevention of gas cooling in low circular velocity haloes 
due to the photoionization of the intergalactic medium, (6) 
the loss of galaxy orbital energy through dynamical friction, 
(7) the subsequent merger between galaxies, which may be 
accompanied by a burst of star formation and (8) chemical 
evolution of the stars and gas. By following these processes, 
GALFORM predicts the star formation history of each galaxy 
(see Baugh 2006 for a review). The code can be run quickly 
for a representative sample of dark matter haloes, making 
it ideal for generating predictions for number counts and to 
populate large volumes to predict galaxy clustering. 

The galaxy formation model we use in this paper is 
that of Baugh et al. (2005; see also Lacey et al. 2010). This 
model reproduces the observed number counts and redshift 
distribution of galaxies at 850 /im, and also the luminosity 
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Table 1. The frequencies covered by the Planck instruments. The rows give: (1) The name of the instrument. The first three frequencies 
(columns 2-4) correspond to the Low Frequency Instrument (LFI) and columns 5-10 to the High Frequency Instrument (HFI). (2) The 
central frequency of the channel in gigahertz. (3) The central wavelength in millimetres. (4) The flux limit for point sources in Janskys 
(taken from Vielva et al. 2003 in the case of LFI and the Planck Collaboration 2011b for the HFI). (5) The angular resolution in 
arcminutes. 
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Frequency (GHz) 
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44 


70 
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143 


217 


353 
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10 


6.81 


4.29 
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2.1 


1.38 
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0.55 


0.35 


Flux limit (Jy) 


0.23 


0.25 


0.24 


0.25 


0.25 


0.16 


0.33 


0.54 


0.71 


Angular resolution (arcmin) 


33 


24 


14 


9.5 


7.1 


5.0 


5.0 


5.0 


5.0 



function of Lyman break galaxies (see Lacey et al. 2011). 
(Note that the more recent model of Bower et al. 2006 does 
not enjoy these successes, and so is not used in this paper.) 
The background cosmology is a spatially flat ACDM model 
(see later for the values of the cosmological parameters) . The 
merger histories of the dark matter haloes are generated 
using an improved Monte Carlo technique that has been 
calibrated against merger trees extracted from an N-body 
simulation (Parkinson, Helly & Cole 2008). 

The model follows two modes of star formation, a "qui- 
escent" mode which takes place in galactic discs and a 
"burst" mode which is triggered by major and minor galaxy 
mergers. Mergers are classified according to the ratio of 
the mass of the merging satellite, M sat , to that of the cen- 
tral galaxy, M ccn . If M sat /M con > / c iii p (where / c ii ip is a 
model parameter) then the merger is defined as a major 
merger. In this case, any cold gas in the two galaxies takes 
part in a starburst which adds stars to the spheroid. Mi- 
nor mergers which have mass ratios in the range /burst < 
Mjat/Mccn < /ciiip and where the primary is also gas rich, 
with M co id/ (M» + M co id) > /gas are also assumed to trigger 
starbursts. In the Baugh et al. model the parameter values 
were set to /cin P = 0.3, /burst = 0.05 and / gas = 0.8. Dif- 
ferent stellar initial mass functions (IMF) are adopted in 
the two modes of star formation. Quiescent star formation 
is assumed to take place with a solar neighbourhood IMF 
(Kcnnicutt 1983). Bursts of star formation are assumed to 
form stars with a top heavy IMF. Baugh et al. (2005) argued 
that the adoption of a top heavy IMF in bursts is necessary 
to reproduce the observed number counts and redshift dis- 
tributions of the faint sub-mm galaxies, whilst at the same 
time reproducing the properties of the local galaxy popula- 
tion. 



2.1.2 The GRASIL spectrophotometry code 

The frequencies sampled by the Planck Low Frequency 
Instrument (LFI; the Planck Collaboration 2011c) and 
High Frequency Instrument (HFI; the Planck Collaboration 
2011d) are listed in Table 1. At these frequencies (corre- 
sponding to wavelengths 0.3-10mm) we assume that the con- 
tribution from galaxies is dominated by star-forming galax- 
ies through dust heated by stars and radio emission resulting 
from gas ionized by stars and synchrotron radiation from rel- 
ativistic electrons accelerated in supernova remnant shock- 
waves. We do not attempt to model the contribution to the 
CBL of active galactic nuclei or galaxies in which the radio 
emission is powered by accretion onto a central black hole. 



To predict the emission from model galaxies at these fre- 
quencies we use the GRASIL spectrophotometric code (Silva 
et al. 1998; Bressan et al. 2002). GRASIL computes the emis- 
sion from the composite stellar population of the galaxy, us- 
ing theoretical models of stellar evolution and stellar atmo- 
spheres. The interaction of the starlight with dust is followed 
with a radiative transfer calculation which assumes a two- 
phase dust medium, and gives the distribution of dust tem- 
peratures within each galaxy using a detailed grain model. 
The output from GRASIL is the galaxy SED from the far-UV 
to the radio (wavelengths 0.01 /xm < A < lm). The main 
features of the hybrid GALFORM-GRASIL model are described 
in Lacey et al. (2010; see also Granato et al. 2000). 



2.2 Populating an N-body simulation with 
galaxies 

We combine the hybrid GALFORM-GRASIL galaxy formation 
model with a high resolution N-body simulation of the clus- 
tering of matter. This allows us to make accurate predic- 
tions for the clustering of galaxies selected by their emission 
at infrared, millimetre and radio wavelengths. In particular 
the small-scale clustering measured in the simulation can be 
significantly different in practice from simple analytical ex- 
pectations based on linear perturbation theory (as shown, 
for example, by Benson et al. 2000). We shall see later that 
for some frequencies, the clustering of galaxy pairs within 
the same dark matter halo is important for the power spec- 
trum of the microwave background intensity fluctuations on 
the angular scales probed by Planck. 

We now describe how galaxies are implanted into an N- 
body simulation. Our starting point is a set of model galax- 
ies generated using the hybrid GALF0RM plus GRASIL code 
set in the concordance ACDM cosmology, which we refer to 
as the MCGAL catalogue. The end point is a galaxy cat- 
alogue implanted in the Millennium Simulation of Springel 
et al. (2005), an N-body simulation of a cosmologically rep- 
resentative volume. We denote this as the MILLGAL cata- 
logue. Some minor complications arise because the cosmol- 
ogy adopted in our original calculation is not quite the same 
as that of the N-body simulation; this issue is dealt with in 
step 2 below. The most accurate calculation of galaxy clus- 
tering is made using the MILLGAL catalogue, so we focus on 
this model in the main paper. In some instances we show 
predictions made using both versions of the model to illus- 
trate regimes in which the MILLGAL predictions are superior. 
The MILLGAL calculation has a finite resolution; as we show 
in the Appendix, this does not affect our results. 
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Figure 1. The halo occupation distribution (HOD) of galaxies 
in the MCGAL catalogue at different redshifts, as indicated by 
the labels. The blue lines show the HOD of central galaxies and 
red lines show the HOD of satellite galaxies. The vertical line 
indicates the halo mass resolution of the Millennium simulation 
and the horizontal line shows (N) = 1. All galaxies with stellar 
mass in excess of 10 7 h^ 1 Mq are allowed to contribute to the 
HOD plotted. This is the input HOD used to build the MILLGAL 
catalogue. 



The MCGAL catalogue is constructed by sampling 
galaxies according to their stellar mass from a much larger 
catalogue of galaxies. This larger catalogue is generated from 
halo merger histories generated using a Monte Carlo method 
(Cole et al. 2000; Parkinson et al. 2008). Star formation 
histories are extracted for the selected galaxies, which are 
then input to GRASIL, along with the predicted dust masses 
and sizes of the disc and bulge components, to compute the 
galaxy SED. For each galaxy we have the stellar mass, host 
halo mass, a weight based on the halo number density and 
on the sampling rate as a function of stellar mass, and the 
galaxy SED. Due to the computational expense of generat- 
ing the MCGAL catalogue, and its similarity with the MILLGAL 
catalogue which we wish to build, we have devised the fol- 
lowing scheme to take advantage of the availability of this 
calculation, rather than making a new calculation based on 
star formation histories extracted directly from galaxies in 
the Millennium simulation. 

The steps followed to populate the Millennium simula- 
tion with galaxies starting from the MCGAL catalogue are 
as follows: 
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Figure 2. The dark matter halo mass function at Z = 0.1. The 
red dotted line shows the halo mass function used in the MCGAL 
catalogue. The black dashed line shows the halo mass function in 
the Millennium simulation at this redshift. The black solid line 
shows the rescaled MCGAL halo mass function, after applying 
a single mass independent adjustment to the halo mass in the 
MCGAL cosmology. 



central (blue line) and satellite (red line) galaxies separately. 
Fig. 1 shows that the HODs of central galaxies in the MC- 
GAL catalogue are approximately step functions. The shape 
of the satellite galaxy HOD is close to a power law in halo 
mass. The lowest mass dark matter halo which appears in 
the MCGAL HODs varies with redshift. The halo mass grid 
used in this calculation is defined at each redshift in order 
to sample haloes with a representative range of abundances. 
At z = 0.1, the lowest mass halo considered is nearly the 
same as the lowest mass halo which can be resolved in the 
Millennium simulation (~ 10 10 ' 3 /i _1 Mq, which is shown by 
the vertical dotted line in Fig. 1). The lowest mass dark 
matter halo in the Millennium simulation is the same at all 
redshifts. We cannot transplant galaxies from the MCGAL 
catalogue which reside in haloes below the mass resolution 
of the Millennium. We test the sensitivity of our predictions 
to this limitation of the N-body catalogue in the Appendix. 
(2) Match the halo mass function between the cosmologies 
used. The MCGAL catalogue, for historical reasons, assumes 
a slightly different cosmology to that adopted in the Mil- 
lennium simulation. The Millennium cosmology is based on 
the first year of WMAP observations 1 . At a given mass, 
the number density of haloes is slightly different in the two 



(1) Construct the halo occupation distribution (HOD) of 
galaxies. GALFORM predicts how many galaxies are contained 
within each dark matter halo. The HOD quantifies the mean 
number of galaxies per halo as a function of the halo mass 
(Benson et al. 2000; Peacock & Smith 2000; Berlind & Wein- 
berg 2002) . Fig. 1 shows the HOD of galaxies in the MCGAL 
catalogue at three different redshifts. We plot the HOD for 



1 The cosmological parameters in the Millennium are: matter 
density Cm = 0.25, cosmological constant = 0.75, Hubble 
constant Ho = 73 km s _1 Mpc~ 1 , primordial scalar spectral in- 
dex n B = 1, baryon density f2t, = 0.045 and fluctuation amplitude 
as = 0.9. The parameters in the MCGAL case are Om = 0.3, 
C A = 0.7, H = 73kms~ 1 Mpc -1 , n s = 1, fi b = 0.04, and 
as = 0.93 (Baugh et al. 2005). 
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cosmologies. The GALFORM-GRASIL code is computationally 
expensive to run, which prohibits running the calculation 
again in the Millennium cosmology. A re-run would also re- 
quire some retuning of the galaxy formation parameters. To 
transplant the galaxies from the halo population in one cos- 
mology to that in the other cosmology, we instead chose 
to relabel the halo masses in the MCGAL model to force a 
match with the Millennium simulation mass function. Fig. 2 
shows the halo mass functions in the original MCGAL calcu- 
lation (red dotted line) and in the N-body simulation (black 
dashed line), along with the rescaled MCGAL halo mass 
function (black solid line) at z = 0.1. To match the halo 
mass functions at z = 0.1 we reduce the halo mass glob- 
ally in the MCGAL catalogue by a factor of 1.2. The halo 
mass function of the rescaled MCGAL catalogue and that 
of Millennium simulation agree to better than 5% over four 
decades in halo mass. We apply the same scheme to other 
redshifts and find that it works equally well, but with slightly 
different scaling factors. The small difference between the 
halo mass functions in the two cosmologies supports our de- 
cision not to rerun the GALFORM-GRASIL calculation in the 
cosmology of the Millennium simulation. 
(3) Place galaxies in the N-body simulation. We generate 
the MILLGAL catalogue of galaxies using the MCGAL 
HOD with halo masses rescaled, as explained in Step 2. The 
central galaxy is placed at the centre of mass of the halo. In 
the case of halo masses for which the HOD specifies N < 1, 
a fraction of the haloes of this mass are populated with a 
single galaxy at random with a probability N. The number 
of satellite galaxies in a given halo mass is assumed to have 
a Poisson distribution for N > 1, with the mean number 
of satellite galaxies given by the HOD. Satellite galaxies are 
assigned to randomly selected dark matter particles which 
are part of the friends-of-friends halo. 



2.3 Calculation of the angular power spectrum of 
intensity fluctuations 

In this subsection we give the theoretical background to 
the calculation of the angular power spectrum of the in- 
tensity fluctuations of galaxies. The discussion is split into 
three parts: §2.3.1 The calculation of the angular correlation 
function of intensity fluctuations. §2.3.2 The calculation of 
the angular power spectrum of intensity fluctuations. §2.3.3 
The estimation of the spatial correlation function of inten- 
sity fluctuations. The first two parts are completely general. 
In the final section we outline how the luminosity-weighted 
spatial correlation function is estimated in the two cases we 
consider: the direct, simulation based approach (MILLGAL) 
and the analytical calculation (MCGAL). Throughout we 
discuss various quantities which depend on intensity at a 
particular frequency. For ease of reading, we suppress the 
explicit frequency dependence in our notation. For exam- 
ple, we write the luminosity density at position x, pi(x, v) 
as Pl(x). We remind the reader that all quantities which 
depend on intensity or luminosity have a frequency depen- 
dence. 



2.3.1 Calculation of the angular correlation function of 
intensity fluctuations 

We can define the spatial correlation function of luminosity 
density, £l(x), as 



(pl(xi)pl(x 2 )) = (pl) 2 (1 + £l(xi - x 2 )) , 



(1) 



where pi(xi) is the luminosity density at position x\ and 
(pl) is the mean luminosity density. Note that here we ne- 
glect shot noise (see §2.3.2) by assuming x\ 7^ x 2 . The cor- 
relation function of galaxy luminosity is related to the stan- 
dard spatial correlation function through 



J f (Li, L 2 , r)L 1 L 2 n(L 1 )n(L 2 )ALiAL 2 
(Jn(L)dL) 2 ' 



(2) 



where n(L) is the mean number density per unit luminosity 
of objects with luminosity L, L t is the luminosity of the i th 
galaxy in the pair, and ^(Li,L 2 ,r) is the cross-correlation 
function of galaxies with luminosities L\ and L 2 . We can 
write the luminosity density as (pi) = f n(L)dL. Similarly, 
we can define an angular surface brightness correlation func- 
tion, wi(6), as 



(/(£)/(&)) = (I) 2 (l + wi(9i - 62)) 



(3) 



where 6 is in radians and the surface brightness is related to 
the comoving luminosity density via 



1(6) d 2 6 



X 2 pL 

dl(x) 



dxd 2 e, 



(4) 



where we have assumed a geometrically flat universe 
(fitotai = 1), and hence di(x), the luminosity distance to 
comoving distance x, is given by di(x) = (1 + z)x. The 
above equation applies in the case of bolometric luminosity 
densities and intensities. However, in practice we are nearly 
always interested in fluxes that are measured over a lim- 
ited frequency band. This introduces an extra (1 + z) factor 
to account for the change in the band width with redshift, 
giving an intensity per unit frequency of 



10) d 2 e = -!- 



{l+ j]i pL dxd 2 e, 

d\(x) 



(5) 



where, if v is the observed frequency, then pi is now the 
luminosity density per unit frequency in a band centred on 
the rest-frame frequency u(l + z). (NB the frequency depen- 
dence of the intensity is suppressed in our notation.) With 
the assumption of a spatially flat universe, the mean inten- 
sity is given by: 



(I) 



^ Jo 



(pl(x)) 
(1 + *) 



da;. 



(6) 



We use Limber's approximation to obtain an expression 
for the angular clustering of flux from the spatial correla- 
tion function (Limber 1953). First, it is assumed that the 
mean number density of galaxies, (n(x)), varies sufficiently 
slowly with redshift (here labelled by the comoving radial 
coordinate x) that over the range of pair separations for 
which £(a?i — x 2 ) 7^ 0, (n(xi)) ~ (n(x 2 )). Second, we assume 
the small angle approximation, i.e. the angular separation 
of pairs of galaxies for which £(a?i — x 2 ) 7^ is small (i.e. 
\6i — 2 \ <C 1 with 6 in radians). 
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Using the above approximations, we can relate the spa- 
tial correlation function of galaxies, £(r), to the angular cor- 
relation, w(8), through Limber's equation 

rz^W)W+*V) 1/2 )d*dn 

w(6) = 2 > (') 

(J" x 2 n(a;)da;) 

where u is a comoving separation parallel to the line of sight, 
such that r 2 = u 2 +x 2 9 2 (again, for small angle separations). 
The analogous relation to Limber's equation for wi(6) is 
given by 



(8) 



{i + zy 



£ L ((u 2 +x 2 6 2 ) 1/2 )dxdu. 



The correlation function of intensity can be evaluated at 
discrete redshifts to give {pl(x)) and £i(r, x) (where redshift 
is again labelled by radial comoving distance x) which can 
then be input into Eq. 8 to compute the angular clustering 
of the flux. Note that in the calculations presented in this 
paper we are interested in the galaxies with flux fainter than 
the detection limits in the Planck bands, as listed in Table 
1. We test the impact of the finite resolution of the N-body 
sample on our predictions in the Appendix. The quantities 
we need to calculate are (7) at the frequencies corresponding 
to the Planck bands and the fluctuations in this background, 
which are given by (I) 2 wi(9). These quantities are predicted 
by the galaxy formation model described in the previous 
subsections. 



2.3.2 Calculation of the angular power spectrum 

The angular power spectrum of the intensity fluctuations 
can be obtained from the angular correlation function of 
intensity using 



{' TV 



Ci{0) = 2tt</) 2 / w/(0)Pi(cos0)sin0d0. 



(9) 



The discrete nature of the sources contributing to the 
CBL, even though they may not be resolved individually 
by an instrument such as Planck, leads to shot noise in the 
intensity fluctuation power spectrum. Even if the galaxies 
displayed no intrinsic clustering and were distributed at ran- 
dom, they would make a contribution to the power spectrum 
through the shot noise. The contribution to the power spec- 
trum from the shot noise depends on the number of sources 
(Tegmark & Efstathiou 1996): 



(10) 



where S hm is the flux detection limit in a given band. 



2.3.3 Estimation of the intensity-weighted correlation 
function 

We follow two different approaches to estimate the spatial 
clustering of galaxies, depending on whether we are using the 
MILLGAL catalogue (for which we have galaxy positions) 
or the MCGAL catalogue (for which we know the mass of 
the halo which hosts each galaxy). In both cases, there is an 



assumption that the clustering of haloes in which galaxies 
are placed depends only on their mass and not on their envi- 
ronment (for an assessment of how halo clustering depends 
on properties besides mass see e.g. Gao et al. 2005; Angulo 
et al. 2008). 

The primary method is to compute the spatial 
luminosity-weighted galaxy correlation function directly, us- 
ing the galaxies transplanted into the N-body simulation 
(i.e. the MILLGAL catalogue). The correlation function is 
estimated from the pair counts of galaxies. This method 
naturally accounts for any difference between the clustering 
of the galaxies and the underlying dark matter because of 
the imposition of the HOD. This approach also allows an ac- 
curate prediction of the correlation function on small scales, 
corresponding to galaxy pairs within the same dark matter 
halo. 

The second approach (used in conjunction with the MC- 
GAL catalogue) is analytical and is intended to show on 
which scales the improvement comes from using the N-body 
simulation. In this case there are two steps in the calculation. 
The first is step is to compute the effective clustering bias 
of the galaxy sample, using the HOD to perform a weighted 
average of the bias of each galaxy based on the analytical 
halo bias (Sheth, Mo & Tormen 2001; see Kim et al. 2009 
for the steps connecting the HOD to the effective bias). The 
analytic estimate of the clustering accounts for only the 2- 
halo term, and assumes that at large separations, the bias 
is independent of scale and depends only on halo mass. The 
effective bias of a galaxy sample is given by 



?L,eff = 



/ / b(M)LN(L\M)n(M)dLdM 
J J LN(L\M)n(M)dLdM ' 



(11) 



where b(M) is the clustering bias factor for haloes of mass 
M, N(L\M) is the HOD (the mean number of galaxies per 
halo which satisfy the sample selection i.e. have luminosity 
L) as a function of halo mass M and n(M) is the halo mass 
function, which gives the abundance per unit volume per 
din M bin of haloes of mass M. 

The second step is to generate a correlation function for 
the dark matter. This is done by Fourier transforming the 
nonlinear power spectrum of the mass derived using the ap- 
proximate analytical prescription of Smith et al. (2003) . The 
galaxy correlation function is then obtained by multiplying 
the nonlinear matter correlation function by the square of 
the effective bias. 



3 RESULTS 

In this section we present results obtained using both the 
MCGAL and MILLGAL catalogues. We also show analytic 
predictions for the clustering of SFGs, and compare these 
to the more accurate predictions obtained from the MILL- 
GAL catalogue, to highlight the shortcomings of the analyt- 
ical calculations. We start by showing the predicted number 
counts of galaxies at the Planck frequencies to show how well 



2 The simulation volume is periodic so the volume of the spherical 
shell for pair separations in the range r to r+dr can be calculated 
analytically; see e.g. Eke et al. 1996 for the form of the estimator 
for the two-point correlation function in this case. 
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Figure 3. The predicted differential number counts from the MCGAL model (solid lines) compared with observational counts (points) 
in the Planck HFI bands, in some cases measured by other telescopes at similar wavelengths. Note that we have adopted the cosmology of 
the MILLGAL model, rescaling the halo masses as described in Sec. 2.2 and using the abundance of haloes in the MILLGAL cosmology. 
The observational data comes from: at 217 GHz, the South Pole Telescope (Vieira et al. 2010; squares) and the Planck Collaboration 
2011e (circles); at 353 GHz from Coppin et al. (2006); and at 545 GHz and 857 GHz from the Herschel- ATLAS Science Demonstration 
Phase field (Clements et al. 2010). The counts are multiplied by flux to the power 2.5 to allow the contrast with the Euclidian counts to 
be better appreciated. 



the model reproduces the observed sub-millimetre counts 
(Sec. 3.1). Next we look at the contribution to the lumi- 
nosity density from galaxies at different redshifts (Sec. 3.2). 
In Section 3.3 we show the predictions for the effective bias 
and contrast the analytical and direct estimates of the spa- 
tial correlation function of intensity. Finally, in Section 3.4 
we present the main predictions of the paper for the angular 
power spectrum of fluctuations in the CBL. 



3.1 The number counts of galaxies at 
sub-millimetre wavelengths 

We first show the model predictions for the number counts of 
galaxies in the Planck HFI bands for which we later present 
clustering predictions. Fig. 3 shows the predicted differential 
galaxy counts and compares these with recent observational 
estimates. The top-right panel of Fig. 3 is an update of the 
comparison shown by Baugh et al. (2005) who considered 
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Figure 4. The luminosity density in the Planck wavebands as a 
function of redshift, predicted using the MCGAL catalogue. 



Figure 5. The contribution to the mean intensity per unit red- 
shift interval in the Planck wavebands in the MCGAL catalogue 
(see Eq. 6). 



the 850 /im number counts. Baugh et al. showed that the 
MCGAL model reproduces the number counts and redshift 
distribution of 850 /im selected samples (see also Almeida 
et al. 2011). The recent measurements of the number counts 
at 250 /im, 350 pm and 500 using the Science Demon- 
stration Phase data from the Herschel Space Telescope by 
Clements et al. (2011) and Oliver et al. (2011) suggest that 
the model predicts too many sources at bright fluxes (see 
predictions in Lacey et al. 2010). This is apparent from the 
comparison between model and observations in the other 
panels of Fig. 3. 

The predictions for fluctuations in the CBL are sensitive 
to the flux-weighted abundance of galaxies. For example, 
the expression for the shot noise contribution to the angular 
power spectrum (Eq. 12) depends on the square of the flux. 
This is similar to the weighting of S 2 ' 5 applied to the differ- 
ential number counts in Fig. 3 (which is standard practice in 
the literature to expand the dynamic range of the counts). 
The dominant contribution to the shot noise is from fainter 
fluxes. The model predictions agree best with the observed 
counts at faint fluxes. The model overpredicts the counts 
at bright fluxes at 350 ^m (857 GHz) and 550 ^m (545 GHz) 
and underpredicts the counts at bright fluxes at 1380 pm 
(217 GHz), which is due to the neglect of radio galaxies in 
the model. 



3.2 The luminosity density of galaxies in the 
Planck bands 

The luminosity density (pL ; Eq. 10) in the Planck frequency 
bands is plotted in Fig. 4 as a function of redshift, com- 
puted using all of the galaxies in the MCGAL catalogue. 
At all nine Planck frequencies, the luminosity density in- 
creases from the present day up to z « 2-5 and then stays 
approximately constant or declines gently to z = 10. At a 



given redshift, the amplitude of the luminosity density in- 
creases with observer frame frequency because L in the rest 
frame increases due to the shape of the SED (until the rest 
frame frequency moves past the peak in the dust emission 
spectrum). For a given observer frame frequency, the lu- 
minosity density is a combination of the SED sampled in 
the rest frame and the abundance of galaxies emitting at a 
given flux. The overall shape of the luminosity density as a 
function of redshift therefore tracks the star formation rate 
density in the universe, with a much gentler decline to high 
redshift, due to the increase in the rest frame uL (due to 
the negative fc-correction) offsetting the overall drop in star 
formation density. 

The fraction of the overall luminosity density that is 
contributed by high redshift sources drops with increas- 
ing frequency, in agreement with previous interpretations 
of fluctuations in the CBL using empirical models (e.g. Hall 
et al. 2010). 

The intensity fluctuation power spectrum depends on 
the mean intensity, which is an integral over the luminosity 
density as given by Eq. 6. From this equation, the contribu- 
tion to the mean intensity per unit redshift interval is given 
by: 

d ffl 1 PL Ax 

dz 4n(l + z)dz' 1 ' 

By plotting this quantity, we can see which redshifts con- 
tribute most to the mean intensity. Fig. 5 shows the con- 
tribution to the mean intensity with redshift in the Planck 
wavebands. A comparison between this plot and Fig. 4 shows 
that the mean intensity is dominated by lower redshifts than 
the luminosity density. 
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Figure 6. The two point luminosity correlation function at 
z = 0.1 (black lines) and z = 2 (red lines) for galaxies fainter 
than the Planck flux limit at 30 GHz. The direct estimates from 
the MILLGAL catalogue are shown by solid lines. The analytic 
calculations derived from the MCGAL sample are shown by dot- 
ted lines when using all galaxies and dashed lines when only using 
those galaxies which reside in haloes above the resolution limit 
of the Millennium. Note that the dotted and dashed lines coin- 
cide with one another at both redshifts, indicating that the finite 
resolution of the Millennium simulation has little impact on the 
predicted clustering of luminosity. 

3.3 The clustering of SFGs in the Planck bands 

We now turn our attention to the computation of the flux 
correlation function. The use of the Millennium N-body sim- 
ulation allows us to make accurate predictions for the clus- 
tering of galaxies on small and intermediate scales. This is il- 
lustrated in Fig. 6 which contrasts the direct estimate of the 
luminosity correlation function estimated from the MILL- 
GAL catalogue with the analytic calculation based on the 
MCGAL catalogue. Recall that the pairwise galaxy counts 
are weighted by luminosity here. At small pair separations, 
r < l/i _1 Mpc, the N-body estimates are up to an order of 
magnitude larger than the analytic ones. Even though the 
analytic calculation takes into account the nonlinear form of 
the matter correlation function, the galaxy correlation func- 
tion can be significantly different on these scales, since the 
analytic calculation ignores the 1-halo term (Benson et al. 
2000). The differences between the two estimates persist to 
intermediate separations of a few megaparsecs, in the tran- 
sition region from the one-halo to two-halo contributions to 
the correlation function. 

In Fig. 7 we take a step closer to the angular power spec- 
trum of intensity fluctuations by plotting the angular flux 
correlation function multiplied by the square of the mean in- 
tensity fluctuations. We also compare the direct estimate of 
the clustering from the N-body simulation (solid lines) with 
the analytic one (dotted lines), for galaxies fainter than the 
Planck detection limits. The amplitude and shape of the two 
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Figure 7. The product of the angular correlation function of 
intensity fluctuations and the square of the mean intensity for 
undetected galaxies in the nine Planck wavebands. Dotted lines 
are for galaxies in MCGAL catalogue using the analytical calcu- 
lations of the bias factor and correlation function. Solid lines are 
for galaxies in MILLGAL catalogue, for which the clustering of 
luminosity is calculated using the Millennium N-body simulation. 

estimates of the correlation function are nearly the same at 
large angular separations. On small angular scales, e.g. 10~ 3 
degrees, the two estimates differ by up to an order of magni- 
tude, due to the more accurate treatment of the 1-halo term 
and the transition between the 1- and 2-halo regimes in the 
N-body calculation (as seen in Fig. 6). 

3.4 The angular power spectrum of fluctuations 
in the CBL from SFGs 

The main results of the paper are shown in Figs. 8 and 
9. Fig. 8 shows the angular power spectrum of the intensity 
fluctuations of undetected galaxies in the Planck wavebands. 
Different components of the model predictions are shown in 
this plot. The long-dashed line shows the intrinsic clustering 
predicted using the N-body simulation, without any contri- 
bution from shot noise. This is to be contrasted with the 
dotted lines, which show the analytic clustering estimate. At 
small angular scales, I > 3000, the N-body estimate exceeds 
the analytic one. Fig. 8 also lets us compare the amplitude 
of the fluctuations from extragalactic sources to the primor- 
dial CMB signal. In the LFI channels, the primordial sig- 
nal dominates on all angular scales. At the HFI frequencies, 
the extragalactic and primordial signals become comparable 
above a particular value of /. The angular power spectrum 
from undetected faint extragalactic sources exceeds the pri- 
mordial power spectrum from I ~ 1000 at 353 GHz, I ~ 100 
at 545 GHz and for all angular scales at 857 GHz. 

Early measurements by Planck have been used to es- 
timate the fluctuations in the CBL from unresolved extra- 
galactic sources (The Planck Collaboration 2011b). As we 
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Figure 8. The angular power spectrum of the intensity fluctuations of undetected galaxies in the nine Planck wavebands. Different 
colours show the frequency bands as listed in the key. The left panel shows the LFI frequencies and the right panel shows the HFI 
wavebands. The dotted lines show the analytic estimates of the intrinsic clustering, which should be compared with the estimates made 
using the Millennium simulation, which are shown by the long-dashed lines. The thin solid lines show the full prediction for the clustering, 
combining the intrinsic clustering with the shot noise derived from the number counts of unresolved sources. The thick lines show the 
primordial CMB power spectra. The vertical arrows indicate the angular resolution of Planck listed in Table 1 and are colour coded by 
frequency. 



have seen in Fig. 8, the primordial signal is orders of mag- 
nitude larger than the extragalactic signal at the lowest fre- 
quencies measured by Planck. Maps at these frequencies can 
be used to estimate the primordial signal expected in the 
higher frequency bands. These authors remove the Galactic 
foreground using observations of emission from neutral hy- 
drogen to clean thermal dust emission. The estimated fluc- 
tuations in the cosmic infrared background due to the un- 
resolved extragalactic sources are shown by the data points 
in Fig. 9. The error bars show the statistical and systematic 
errors (the area used covers 140 square degrees in six differ- 
ent fields). The systematic error, e.g. due to uncertainties in 
the beam, exceeds the statistical error at high I in the higher 



frequency channels (see the Planck Collaboration 2011b for 
further details). 

The theoretical predictions for the extragalactic inten- 
sity fluctuations in the CBL are considered more closely in 
Fig. 9. The shot noise computed from the MILLGAL model 
using Eq. 10 is shown by the short dashed line. The shot 
noise makes a fixed contribution to Ci but has a scale de- 
pendence in Fig. 9 since here we plot 1(1 + 1)C;. The contri- 
bution to Ci from the intrinsic clustering of the unresolved 
galaxies, as estimated using the N-body based MILLGAL 
catalogue, is shown by the long dashed line. The correspond- 
ing analytical estimate is shown by the dotted line. At the 
smallest scales plotted, the N-body estimate exceeds the less 
accurate analytical one by a factor of three or more. The 
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Figure 9. The angular power spectrum of intensity fluctuations in the extragalactic cosmic infrared and millimetre background. Each 
panel corresponds to a different HFI frequency as indicated. The points show the extragalactic background fluctuation estimated from 
Planck measurements (The Planck Collaboration 2011b), after removing the cosmological signal and the Galactic foreground. The error 
bars show the statistical and systematic errors. The lines show our theoretical predictions. The short-dashed line shows the shot noise 
from the counts predicted by the MILLGAL model. The long dashed line shows the clustering estimated from the MILLGAL catalogue, 
without any shot noise, which we refer to as the intrinsic clustering. The analytic estimate of the intrinsic clustering is shown by the 
dotted line. The solid line shows our full prediction for the intensity fluctuations, combining the intrinsic clustering (long dashed line) 
and the shot noise (short dashed line). Another version of this prediction, derived by combining the intrinsic clustering with the shot 
noise estimated using the model of Bcthcrmin et al. (2011) is shown by the dot-dashed line. 



overall prediction for the intensity fluctuation power spec- 
trum, combining the intrinsic clustering with the shot noise, 
is shown by the solid line. The shot noise makes an im- 
portant contribution to the power spectrum at high I. We 
have commented already that our model overpredicts the 
number counts of galaxies in the Herschel bands at bright 



fluxes. This will result in turn in an overprediction of the 
shot noise in these bands. To illustrate how this can affect 
our predictions, we also show a version of the power spec- 
trum in which we replace the shot noise calculated using our 
model with that from the empirical number count model of 
Bethermin et al. (2011). This alternative prediction is shown 
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by the dot-dashed line in Fig. 9. We note that these predic- 
tions are not meant to supersede those shown by the solid 
line, as now the number of sources is decoupled from the 
calculation of the intrinsic clustering. 

Fig. 9 shows that the model predictions are generally 
within a factor of three or better of the observational es- 
timate of the extragalactic background fluctuations. It is 
notable that in some of the HFI bands, the predicted shape 
of the power spectrum is similar to the observational esti- 
mate. The agreement between the model and observations 
is best at small angular scales in the two highest frequency 
channels. There is a mismatch in amplitude on larger scales 
(smaller I) in these higher frequency channels. 

3.5 How can the model predictions be improved? 

The goal of this paper is to present a new framework for pre- 
dicting the contribution of star-forming galaxies to intensity 
fluctuations in the cosmic background light. We have used a 
previously published model to illustrate the technique. Un- 
like other studies in the literature, we have not varied any 
model parameters in order to improve the agreement of the 
model predictions with the signal inferred from observations. 

Fig. 9 shows the model underpredicts the observed clus- 
tering on large angular scales in the two highest frequency 
HFI channels, whilst giving a reasonable match to the power 
spectrum on small scales. To improve the model predictions 
at 545 GHz and 857 GHz, we would need to increase the ef- 
fective bias of the unresolved sources, whilst reducing their 
number slightly. This requires an increase in the typical ef- 
fective host halo mass. Such a shift could be achieved by 
reducing the efficiency of star formation in galaxies hosted 
by low mass haloes. 

On the other hand in the two lowest frequency chan- 
nels shown in Fig. 9, 217 GHz and 353 GHz, the model over- 
predicts the clustering, so there needs to be a reduction in 
the effective host halo mass in these cases. We note that 
the model works best overall at 353 GHz, the frequency at 
which the model was originally tuned to match the observed 
galaxy number counts and redshift distribution. 

Fig. 10 gives some insight into which haloes dominate 
the clustering signal. The top panel of Fig. 10 shows the total 
star formation rate summing over galaxies in the same dark 
matter halo, plotted as a function of host halo mass. The y- 
axis is on a linear scale so that we can gain an impression of 
which haloes make the most important contribution to the 
overall star formation rate density (we would also need to 
take into account the halo mass function to connect this plot 
to the luminosity density) . The distribution is dominated by 
central galaxies hosted by dark matter haloes with masses in 
the range 10 11 -10 12 Mq. In this paper we study correla- 
tions in intensity, so an important quantity to consider is the 
total luminosity per halo, summing over galaxies within the 
same dark matter halo, which is plotted as a function of host 
halo mass in the bottom panel of Fig. 10. The far-infrared 
luminosity depends on a galaxy's star formation rate, along 
with its dust content and the level of extinction. The dis- 
tribution of luminosity per halo has a similar form to that 
of the star formation rate, with a more pronounced tail to 
higher halo masses. This plot shows that the mass resolution 
of the Millennium N-body simulation is sufficient for mod- 
elling the intensity fluctuations of unresolved galaxies. To 
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Figure 10. Top panel: the total star formation rate summing 
over galaxies within common dark matter haloes as a function of 
host halo mass. Different redshifts are shown by different coloured 
lines, as indicated by the key. The solid lines show the predictions 
for all galaxies and the dashed lines for satellite galaxies only. 
Bottom panel: the total luminosity per halo in the observer frame 
at 353 GHz as a function of host halo mass. Lines styles have the 
same meaning as in the top panel. 



change the clustering predictions of the model, we need to 
move the location of the peak in the distribution plotted in 
Fig. 10, which means finding a way to put central galaxies of 
a given luminosity into different mass haloes, depending on 
the sense of the change required in the two-halo clustering 
term. The tail of satellite galaxies apparent in higher mass 
haloes influences the one-halo clustering term (i.e. pairs of 
galaxies within the same halo). 
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We have carried out a preliminary exploration of param- 
eter space, changing one aspect of the model at a time, with- 
out trying to retain the previous successes of the model. We 
altered the strength of feedback from supernovae and experi- 
mented with deleting satellite galaxies assigned to subhaloes 
which can no longer be resolved in the Millennium Simula- 
tion. These changes produce a similar shift in the model 
predictions at each frequency and so do not produce the 
differential changes that we need to improve the match to 
the observations. A full, multi-dimensional parameter search 
to find a model with an improved match to the intensity 
fluctuations, which also reproduces the previous successful 
matches to observations, is beyond the scope of the current 
paper. 



4 CONCLUSIONS 

Fluctuations in the intensity of the microwave background 
radiation arise from a number of sources: ripples in the den- 
sity of matter in the primordial Universe, emission from the 
interstellar medium in our galaxy, and extragalactic sources, 
such as star- forming galaxies (SFGs), radio galaxies or the 
hot plasma in galaxy clusters. Identifiable sources can be 
removed from intensity fluctuation maps. Sources below the 
confusion limits of the measurements cannot be explicitly 
excised. Their contribution can be removed statistically in 
analyses of the primordial signal or can be isolated to study 
the history of star formation in the Universe and its connec- 
tion to structure formation in the dark matter. 

Here we have introduced a hybrid model which com- 
bines a physical model of galaxy formation with a N-body 
simulation of the clustering of dark matter to predict the 
contribution of SFGs to the intensity fluctuations in the 
CBL. This is the first time that it has been possible to com- 
pute the abundance and clustering of SFGs together in a 
physical model. The model predicts the radio emission from 
star-forming galaxies and the emission from dust heated by 
stars, in a self consistent manner, with the dust grains in 
thermal equilibrium. The amount of heating depends on the 
star formation and chemical enrichment history of the galax- 
ies, and on their dust mass and linear size; all of these prop- 
erties are predicted by the model. Our calculation tells us 
how many galaxies above a stated flux limit should reside 
within dark matter haloes of a given mass. This information, 
in combination with an N-body simulation of the clustering 
of dark matter, allows a direct calculation of the clustering of 
SFGs. In previous calculations based on dark matter haloes, 
the halo occupation distribution was adjusted by hand, with 
no connection between the properties of the host halo and 
the spectral energy distribution of the galaxy. 

We have illustrated this new framework using a pub- 
lished model of galaxy formation (Baugh et al. 2005). This 
model matches the observed galaxy number counts at 353 
GHz (850 pirn) , along with other observations of the galaxy 
population at low and high redshift. We have taken the pre- 
dictions of this model without any adjustments. Hence, the 
example calculation we present is in effect parameter-free, 
since we make no further adjustment to the values of the pa- 
rameters which specify the galaxy formation model. This is 
an important distinction of our work from halo occupation 
distribution modelling of the CBL fluctuations, in which case 



the model parameters are adjusted to improve the match to 
the measured clustering. In view of this, it is impressive that 
our predictions only disagree with the fluctuations in the 
CBL inferred from early Planck data by at worst a factor of 
three. It is also important to bear in mind that the observa- 
tional estimate of the CBL fluctuations is heavily processed, 
and significant contributions from other sources have to be 
removed to isolate the extragalactic signal. Empirical calcu- 
lations in the literature have a more limited scope than our 
model. Such calculations do contain parameters which are 
tuned to improve directly the agreement between the pre- 
dicted and observed clustering. However, in general these 
models do not connect the emission from galaxies to their 
host dark matter haloes (for an exception see e.g. Shang 
et al. 2012). 

Our model agrees best the inferred CBL fluctuations 
at high frequencies on small scales (high I), for which the 
shot noise from discrete sources and the clustering of galax- 
ies within common dark matter haloes dominate, and where 
the isolation of the extragalactic signal is most secure. By 
using an N-body simulation, we are able to make accurate 
predictions for the one-halo as well as two halo contribution 
to the intensity fluctuations. These predictions are signifi- 
cantly different from simple analytical calculations on small 
angular scales. In general the model does less well on larger 
scales, around / ~ 100. This implies that the example model 
predicts the wrong effective bias or two-halo clustering term. 
At 217GHz, the predicted bias is too high by a factor of 1.6. 
At 857GHz, the effective bias predicted is too small by a sim- 
ilar factor (1.7). As the intrinsic clustering dominates over 
the shot noise on these scales, this suggests that some redis- 
tribution of galaxies between dark matter haloes is required 
in the model. This is more difficult to realise than it may 
sound, as these adjustments would have to be made without 
changing the small-scale clustering by much. Of course, if the 
number of sources is also changed, then this will alter the 
shot noise, which mainly affects the amplitude of the small 
scale clustering. Another possible explanation for the dis- 
crepancy between the model predictions and the Planck re- 
sults is that the emissivity of the dust assumed in the model 
may be incorrect. Baugh et al. modified the dust emissivity 
in bursts from the standard value of e tx v~ 2 used in quies- 
cent star formation to e tx v~ 1,5 . This boosts the emission 
at longer wavelengths, and may in part be responsible for 
the excess counts at low frequencies. 

As new observations of SFGs become available through, 
for example, the Herschel Space Observatory, new con- 
straints will be placed on the galaxy formation model which 
underpins our method (Lacey et al. 2008; 2010). Along with 
improved treatments of key model ingredients, such as star 
formation (Lagos et al. 2010), this will allow us to devise new 
models which better match the new observations of SFGs. 
The clustering of intensity fluctuations in the CBL will offer 
an important additional observational constraint that such 
galaxy formation models will be able to exploit. 
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APPENDIX 

The halo mass resolution of the MCGAL catalogue is essen- 
tially arbitrary, as the full memory of the computer can be 
devoted to a single dark matter halo merger history. Further- 
more, the mass resolution can be adjusted with redshift, to 
ensure that a representative sample of halo masses is mod- 



The contribution of star-forming galaxies to fluctuations in the cosmic background light. 15 
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Figure 11. The luminosity function (LF) at 857 GHz (350 jim) at 
z = 0.1 (left) and z = 5 (right). Solid lines shows the full MCGAL 
catalogue and dotted lines show the result from this catalogue 
when restricted to haloes with masses above the resolution limit 
of the Millennium simulation (M halo = 1.72 X 10 10 h _1 Me). Note 
that the dotted and solid lines are coincident in the left panel. The 
dashed line shows the predictions from the MILLGAL catalogue. 
The blue lines show the contribution to the LFs from central 
galaxies. The red lines show the LF of satellite galaxies. 



elled at each epoch. The mass resolution of the MILLGAL 
sample is set by the Millennium N-body simulation and 
is fixed at all redshifts. In this Appendix we compare re- 
sults from the two calculations to demonstrate the impact 
that the finite resolution of the MILLGAL sample has on our 
model predictions. We conclude that the predictions for the 
intensity fluctuations are insensitive to the resolution of the 
Millennium Simulation. 

We first consider the luminosity functions predicted in 
the two catalogues in Fig. 11. At low redshift, the MILL- 
GAL and MCGAL catalogues are in very good agreement 
with one another, both for central and satellite galaxies (i.e. 
at z — 0.1 in the left panel of Fig. 11). The luminosity 
function of the two catalogues differ at high redshift be- 
cause the MILLGAL does not include galaxies hosted by 
haloes below the mass resolution of the Millennium simula- 
tion. However, the MILLGAL catalogue reproduces well the 
luminosity function of the MCGAL catalogue at higher lu- 
minosities, at which the galaxies tend to be hosted by more 
massive dark matter haloes which are resolved in Millennium 
simulation (see the right panel in Fig. 11). 

In Fig. 12, the luminosity density of galaxies hosted by 
dark matter haloes resolved by the Millennium simulation 
is lower than that of all the galaxies in the MCGAL cata- 
logue for z > 4. Galaxies in low mass dark matter haloes 
contribute significantly to the luminosity density at high z. 

In Fig. 13, we plot the integrated luminosity weighted 
bias. The result shows that there is little impact on this mea- 
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Figure 12. The luminosity density in the Planck wavebands as a 
function of redshift, predicted using the MCGAL catalogue. The 
solid lines show the predictions using all galaxies in the MCGAL 
catalogue and the dotted lines show the results using only those 
galaxies hosted by haloes which could be resolved in the Millen- 
nium Simulation. 

sure of the intensity fluctuations on applying the resolution 
limit of the N-body simulation to the MCGAL catalogue. 

The similarity between the results shown in Fig. 13 for 
different halo mass resolution limits and the fact that the 
mean intensity is dominated by low redshifts imply that our 
predictions for the power spectrum of intensity fluctuations 
should be insensitive to the resolution limit of the MCGAL 
catalogue. This is confirmed in Fig. 14 in which we plot 
the analytic estimate of the product of the square of the 
mean intensity and the angular correlation function of in- 
tensity using the MCGAL catalogue. The predictions for 
the full MCGAL catalogue are indistinguishable from those 
restricted to galaxies which could be resolved in the Millen- 
nium simulation. 
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Figure 13. The integrated luminosity weighted bias (the numer- 
ator of the effective bias as defined in Eq. 11) as a function of 
redshift in the Planck wavebands calculated using the analytic 
halo bias from Sheth et al. (2001). Only galaxies fainter than the 
Planck detection limits listed in Table 1 contribute. The solid 
curves show the predictions using all the galaxies in the MCGAL 
catalogue. The dotted curves show the results using only those 
galaxies hosted by haloes which could be resolved in the Millen- 
nium simulation. 
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Figure 14. The product of the angular correlation function of 
intensity fluctuations and the square of the mean intensity for 
undetected galaxies in the nine Planck wavebands. Solid lines 
are for all galaxies in MCGAL catalogue and dotted lines are 
for galaxies in haloes which could be resolved in the Millennium 
simulation, using the analytical calculations of the bias factor and 
correlation function. These predictions are very similar, showing 
that halo mass resolution has little impact on our predictions. 



